! The following program creates clinoform surfaces (see Chapter 2 section 2.3)
! and corresponding heterogeneity along those surfaces (see Chapter 3 section 3.7.3)
! Here we briefly desribe the subroutines used throughout this program
!----------------------------------------------------------------------------------------------------------------------
! The first section of algorithm initiates the arrays that are used in subsequent subroutines
program clinoform_modelling

   implicit none
   
   interface
    subroutine bounding_surfaces_read_in(BOUNDING, n, m, xmin, xmax, ymin, ymax, nonvalue, base_name, top_name)
    character :: base_name*1000, top_name*1000
    integer :: n, m
    real :: xmin, xmax, ymin, ymax, nonvalue
    real, allocatable :: BOUNDING(:,:,:)
    end subroutine
   end interface
   
   interface
    subroutine array_allocation (temp_name,TEMP,n,m, nonvalue,xmin, xmax, ymin, ymax)
    character :: temp_name*1000 
    integer :: n,m
    real :: xmin, xmax, ymin, ymax, nonvalue
    real, allocatable :: TEMP(:,:,:)
    end subroutine
   end interface
   
   interface 
    subroutine write_zmap_file(TEMP,temp_name,n,m,k,xmin,xmax,ymin,ymax)
    character :: temp_name*1000 
    integer :: n,m,k
    real :: xmin, xmax, ymin, ymax
    real, allocatable :: TEMP(:,:,:)
    end subroutine
   end interface
    
    interface
     subroutine increase_grid_area(n,m,BOUNDING, nonvalue, xmin, xmax, ymin, ymax)
     integer :: n,m
     real :: xmin, xmax, ymin, ymax, nonvalue
     real, allocatable :: BOUNDING(:,:,:)
     end subroutine
    end interface
    
    interface
     subroutine increase_nodes(n, m, BOUNDING, nonvalue, n_regrid, m_regrid, NODE_CHANGE, numsurf, MORE_NODES_BOUNDING, NODE_CHANGE_1)
     integer :: n_regrid,m_regrid, n, m, numsurf
     real :: nonvalue
     real, allocatable :: BOUNDING(:,:,:)
     real, allocatable :: MORE_NODES_BOUNDING(:,:)
     real, allocatable :: NODE_CHANGE(:,:)
     real, allocatable :: NODE_CHANGE_1(:,:)
     end subroutine
    end interface
    
    interface
     subroutine rotate_surfaces(n_regrid, m_regrid, NODE_CHANGE, numsurf, ROTATION, xmin, xmax, ymin, ymax, user_angle, nonvalue, resolx)
     integer :: n_regrid, m_regrid, numsurf
     real :: nonvalue, xmin, xmax, ymin, ymax, user_angle, resolx
     real, allocatable :: ROTATION(:,:,:)
     real, allocatable :: NODE_CHANGE(:,:)
     end subroutine
    end interface
    
    interface
     subroutine create_clinoform_surfaces(ROTATION, nonvalue, n_regrid, m_regrid, resolx, xmin, xmax, ymin, ymax, user_angle, A, filepath, proportion, Petrel)
     character :: filepath*1000
     integer :: n_regrid, m_regrid, A, Petrel
     real :: nonvalue, resolx, xmin, xmax, ymin, ymax, user_angle, proportion
     real, allocatable :: ROTATION(:,:,:)
     end subroutine
    end interface
    
    interface
     subroutine post_processing(A, processing, base_name, top_name, model, filepath, Petrel)
     character :: base_name*1000, top_name*1000, filepath*1000
     integer :: A, processing, model, Petrel
     end subroutine
    end interface
      
    character :: temp_name*1000,base_name*1000, top_name*1000, filepath*1000
    integer :: n, m, numsurf, n_regrid, m_regrid, k, SPACE, A, processing, model, Petrel
    parameter (SPACE=6100)
    real :: nonvalue, xmin, xmax, ymin, ymax, user_angle, resolx, proportion
    real, allocatable :: BOUNDING(:,:,:)
    real, allocatable :: MORE_NODES_BOUNDING(:,:)     
    real, allocatable :: NODE_CHANGE_1(:,:)
    real, allocatable :: NODE_CHANGE(:,:)
    real, allocatable :: ROTATION(:,:,:)
    
!--------------------------------------------------------------------------------------------------------------------------
!initiates new seed for stochastic modelling of heterogeneity along clinoform surfaces
        call random_seed()
!--------------------------------------------------------------------------------------------------------------------------
!subroutine that creates arrays containing data on the bounding surfaces that define the volume within which clinoforms 
!are modelled (Chapter 2, section 2.3.1)from user specified input surfaces
        call bounding_surfaces_read_in(BOUNDING, n, m, xmin, xmax, ymin, ymax, nonvalue, base_name, top_name)
!--------------------------------------------------------------------------------------------------------------------------
!subroutine that increases the area of the grid on which the input surface are modelled, to allow for rotation of the surfaces 
!in a later triangulation step
        call increase_grid_area(n, m, BOUNDING, nonvalue, xmin, xmax, ymin, ymax)
!--------------------------------------------------------------------------------------------------------------------------
!user specified progradtion direction       
        print*, "Enter the angle clinoform progradation angle relative to North"
        read*, user_angle
        
        if (0.0<=user_angle<=90.0)then;
        user_angle=(90.0 - user_angle)
        else if (90.0<user_angle<=180.0)then;
        user_angle=(90-(user_angle-90.0))+270 
        else if (180.0<user_angle<=270.0)then;
        user_angle=(90-(user_angle-180.0))+180.0
        else;
        user_angle=(90-(user_angle-270.0))+90.0
        end if
        
        write(*,*)"Enter the proportion of the clinoform surface that acts as a barrier to flow"
        read(*,*) proportion
        
        allocate (MORE_NODES_BOUNDING(SPACE,SPACE))
        allocate (NODE_CHANGE_1(SPACE,SPACE))
        allocate (NODE_CHANGE(SPACE,SPACE))
        allocate (ROTATION(2,SPACE,SPACE))
!--------------------------------------------------------------------------------------------------------------------------
!subroutines that increase the resolution of the grid on which the input surface are modelled
!and rotate the user defined bounfding surfaces to match the user defined
!progradation direction of clinoform surfaces (Chapter 2, section 2.3.4)  
  do numsurf=1,2
        call increase_nodes(n, m, BOUNDING, nonvalue, n_regrid, m_regrid, NODE_CHANGE, numsurf, MORE_NODES_BOUNDING, NODE_CHANGE_1)
        call rotate_surfaces(n_regrid, m_regrid, NODE_CHANGE, numsurf, ROTATION, xmin, xmax, ymin, ymax, user_angle, nonvalue, resolx)
  end do
             
         deallocate (MORE_NODES_BOUNDING)
         deallocate (NODE_CHANGE_1)
         deallocate (NODE_CHANGE)
         deallocate (BOUNDING)     
        
        call listing(SPACE)  
!--------------------------------------------------------------------------------------------------------------------------
!subroutine that creates the clinoform surface (contains Eq. 2.1-2.8 from Chapter 2 sections 2.3.2 and 2.3.3)     
        call create_clinoform_surfaces(ROTATION, nonvalue, n_regrid, m_regrid, resolx, xmin, xmax, ymin, ymax, user_angle, A, filepath, proportion, Petrel)
        
        print*, "Do surfaces need post-depositional modification?"
        read*, processing
        
        !print*, "Do surfaces need model extent post-processing?"
        !read*, model
        model = 0
        
        call post_processing(A, processing, base_name, top_name, model, filepath, Petrel)
                             
        open(unit=3, file="listing.txt")  
        close(unit=3, STATUS='DELETE')
        open(unit=30,file="passageone.zmap+")
        close(unit=30, STATUS='DELETE')
        open (unit=60, file="passageone_1.zmap+")
        close(unit=60, STATUS='DELETE')
        open(unit=4, file="Clinoform List.txt")  
        close(unit=4, STATUS='DELETE')
        open(unit=5, file="bounding1.zmap+")  
        close(unit=5, STATUS='DELETE')
        open(unit=6, file="bounding2.zmap+")  
        close(unit=6, STATUS='DELETE')
        
    write(*,*) 'Press return to exit program'
    read(*,*)
    
stop
end program clinoform_modelling
!----------------------------------------------------------------------------------------------------------------------------
!subroutine to read in user defined bouding surfaces
subroutine bounding_surfaces_read_in(BOUNDING, n, m, xmin, xmax, ymin, ymax, nonvalue, base_name, top_name)

  interface
    subroutine array_allocation (temp_name,TEMP,n,m, nonvalue,xmin, xmax, ymin, ymax, passage)
    character :: temp_name*1000 
    integer :: n,m
    real :: xmin, xmax, ymin, ymax, nonvalue
    real, allocatable :: TEMP(:,:,:)
    real, allocatable :: passage(:,:)
    end subroutine
   end interface
   
   interface 
     subroutine write_zmap_file(TEMP,temp_name,n,m,k,xmin,xmax,ymin,ymax)
     character :: temp_name*1000 
     integer :: n,m,k
     real :: xmin, xmax, ymin, ymax
     real, allocatable :: TEMP(:,:,:)
     end subroutine
   end interface 
    
    character :: temp_name*1000,top_name*1000,base_name*1000 
    integer :: k, n, m, SPACE
    parameter (SPACE=6100)
    real :: nonvalue, xmin, xmax, ymin, ymax   
    real,allocatable :: BOUNDING(:,:,:)
    real, allocatable :: TEMP(:,:,:)
    real, allocatable :: passage(:,:)
    
    allocate (passage(SPACE,5))
    allocate (TEMP(2,SPACE,SPACE))
    
    write(*,*)"Enter the file path (including file name) of the base bounding surface:"
    read(*,*) base_name
    temp_name=base_name
    call array_allocation (temp_name,TEMP,n,m, nonvalue,xmin, xmax, ymin, ymax, passage)  
    allocate (BOUNDING(2,(SPACE),(SPACE)))
    k=1
        do i=1,n
            do j=1,m
               BOUNDING(k,i,j)=TEMP(1,i,j)
               if (BOUNDING(k,i,j)>nonvalue.and.BOUNDING(k,i,j)<0) then;
                    BOUNDING(k,i,j)= BOUNDING(k,i,j)*-1
               end if
            end do
        end do
     deallocate (TEMP)

    allocate (TEMP(2,SPACE,SPACE))
    write(*,*)"Enter the file path (including file name) of the top bounding surface:"
    read(*,*) top_name
    temp_name=top_name
    call array_allocation(temp_name,TEMP,n,m, nonvalue,xmin, xmax, ymin, ymax, passage)
    k=2
        do i=1,n
            do j=1,m
               BOUNDING(k,i,j)=TEMP(1,i,j)
                if (BOUNDING(k,i,j)>nonvalue.and.BOUNDING(k,i,j)<0) then;
                    BOUNDING(k,i,j)= BOUNDING(k,i,j)*-1
                end if
            end do
        end do
    deallocate (TEMP)
     
     deallocate (passage)  
     
return
end subroutine bounding_surfaces_read_in
!-------------------------------------------------------------------------------------------------------------------------------
!subroutine that creates arrays from user defined surfaces
subroutine array_allocation (temp_name,TEMP,n,m, nonvalue,xmin, xmax, ymin, ymax, passage)

    character :: temp_name*1000
    character :: dummy1*10, dummy2*10, dummy3*10, dummy4*10
    character :: dummy5*10, dummy6*10, dummy7*10, dummy8*10
    integer :: m, n, lines, rest, i, j, h, raw, mode, k, SPACE
    parameter (SPACE=6100)
    real :: numb, nonvalue, xmin, xmax, ymin, ymax, depth_min, depth_max
    real :: passagefive(5), passagefour(4), passagethree(3), passagetwo(2)
    real :: passageone
    real, allocatable :: passage(:,:)
    real, allocatable :: TEMP(:,:,:)  

   open(unit=1, file=temp_name)
   	 read(1,*) dummy1
	 read(1,*) dummy2
	 read(1,*) dummy3
	 read(1,*) dummy4
	 read(1,*) dummy5
	 read(1,*) dummy6
	 read(1,*) numb, nonvalue
	 read(1,*) m, n, xmin, xmax, ymin, ymax
	 read(1,*) dummy7
     read(1,*) dummy8
     
     do i=1,SPACE
       do j=1,SPACE
          TEMP(1,i,j)=nonvalue
          TEMP(2,i,j)=nonvalue
       end do
    end do
      
     lines=m/5+1
     rest=m-(lines-1)*5
	if (rest==0) then;
    	lines=lines-1
    end if
    
	 do i=1,n
    	h=1
        do j=1,(lines-1)
           read(1,*) passagefive
           passage(h,1:5)=passagefive
           h=h+1
    	end do

		if (rest==1) then;
          read (1,*) passageone
          passage(h,1)=passageone
        else if (rest==2) then;
          read (1,*) passagetwo
          passage(h,1:2)=passagetwo
        else if (rest==3) then;
          read (1,*) passagethree
          passage(h,1:3)=passagethree
        else if (rest==4) then;
          read (1,*) passagefour
          passage(h,1:4)=passagefour
        else;
          read (1,*) passagefive
          passage(h,1:5)=passagefive
		end if
        
		do k=1,m
          raw=k/5
          mode=k-raw*5
          if (mode==0) then;
            mode=5
            raw=raw-1
          end if
          TEMP(1,i,k)=passage(raw+1,mode)          
         end do
     end do
    close(unit=1) 
              
return
end subroutine array_allocation
!-------------------------------------------------------------------------------------------------------------------------------
!subroutine that exports all surfaces from the algorithm
subroutine write_zmap_file(TEMP,temp_name,n,m,k,xmin,xmax,ymin,ymax)

    character :: temp_name*1000
    character :: path*97
    integer :: n, m, i, k
    real :: xmin, xmax, ymin, ymax
    real, allocatable :: TEMP(:,:,:)
          
     open(unit=1,file=temp_name)
        write(1,*)'!'
        write(1,*)'!     ZIMS FILE NAME :  Clinoform_Surface.txt'
        write(1,*)'!     FORMATTED FILE CREATION DATE: xx/xx/xxxx'
        write(1,*)'!     FORMATTED FILE CREATION TIME: xx:xx'
        write(1,*)'!'
        write(1,1)'@Grid HEADER, GRID, 5'
1       format(A21)
        write(1,*)'	15,-9999.000,  ,	4 ,1.0000'
        write(1,"(I6,A,I5,A,1F12.4,A,1X,1F12.4,A,1F12.4,A,1F12.4)") m,',',n,',',xmin,',',xmax,',',ymin,',',ymax
        write(1,*)'0.0000     0.0000     0.0000'
	    write(1,2)'@'
2       format(A1)
        do i = 1, n 
            write(1,"(1F10.4,1X,1F10.4,1X,1F10.4,1X,1F10.4,1X,1F10.4)") TEMP(k,i,1:m)
         end do
    close(unit=1)    

return
end subroutine write_zmap_file
!-------------------------------------------------------------------------------------------------------------------------------
!corresponding subroutine that increases the area of the grid on which the input surface are modelled, to allow for rotation of the surfaces 
!in a later triangulation step
subroutine increase_grid_area(n,m,BOUNDING, nonvalue, xmin, xmax, ymin, ymax)

    character :: temp_name*1000
    integer :: SPACE, n, m, icoordinate, jcoordinate, c, i, j, short_side, long_side, resol_y_int, resol_x_int, diagonal_side
    integer :: nodes_req, y_length_int, nodes_actual, nodes_extra, m_original, m_temp_start, m_temp_end, m_mid, n_original
    integer :: n_temp_start, n_temp_end, n_mid, extent, extent_req, k
    parameter (SPACE=3100)
    real :: nonvalue, xmin, xmax, ymin, ymax, icornerC, jcornerC, icornerD, jcornerD, icornerB, jcornerB, last_value    
    real, allocatable :: BOUNDING(:,:,:)
    real, allocatable :: TEMP_2(:,:,:) !4, space, space
       
    allocate (TEMP_2(4,SPACE,SPACE))
!first find first point in bounding surfaces array, xmin
     icoordinate =0
     jcoordinate =0
     c=1
  do while(icoordinate<1)
     i=(c/(m))+1
     j= c-((i-1)*(m-1))

    if (j<1)then;
     j=n
    end if

    if(BOUNDING(1,i,j)>(nonvalue+1))then;
     icoordinate=i
     jcoordinate=j
    end if
     c=c+1
  end do
 
     icornerC = icoordinate
     jcornerC = jcoordinate


!need to calculate last value in file 
  do i=1, n
     do j=1,m
        if(BOUNDING(1,i,j)>(nonvalue+1))then;
          last_value=BOUNDING(1,i,j)
          icornerD=i
          jcornerD=j 
        end if
     end do
  end do
  

     icoordinate =0
     jcoordinate =0
     c=1
  do while(icoordinate<1)
     j=(c/(n-1))+1
     i= c-((j-1)*(n-1))

    if (i<1)then;
     i=n
    end if

    if(BOUNDING(1,i,j)>(nonvalue+1))then;
     icoordinate=i
     jcoordinate=j
    end if
     c = c+1
  end do
  
     icornerB =icoordinate
     jcornerB = jcoordinate
     
!Calculate how many nodes are actually needed
     short_side=SQRT(((icornerB-icornerC)**2)+((jcornerB-jcornerC)**2))+1
     long_side= SQRT(((icornerD-icornerB)**2)+((jcornerB-jcornerD)**2))+1
     
     resol_y_int = ((ymax-ymin)/m)+1
     resol_x_int = ((xmax-xmin)/n)+1
     short_side=short_side*resol_y_int
     long_side=long_side*resol_y_int
     
     diagonal_side=((long_side**2)+(short_side**2))+1
     diagonal_side=diagonal_side**0.5  
     nodes_req=diagonal_side/resol_y_int
    
     if (MOD(nodes_req,2)>0)then;
     nodes_req=nodes_req+1
     end if

     y_length_int=(ymax-ymin)
     nodes_actual=((ymax-ymin)/resol_y_int)+1
     
     
        !increase nodes section
 if (diagonal_side>y_length_int)then;
     
        nodes_extra=nodes_req-nodes_actual
     
        if (MOD(nodes_extra,2)>0)then;
        nodes_extra=nodes_extra+1
        end if  
     
       ymax=ymax+(resol_y_int*nodes_extra)
       ymin=ymin-(resol_y_int*nodes_extra)
       xmax=xmax+(resol_x_int*nodes_extra)
       xmin=xmin-(resol_x_int*nodes_extra)
         
       m_temp_start=nodes_extra
       m_mid=m+nodes_extra
       m_temp_end=m+(nodes_extra*2)
   
       n_temp_start=nodes_extra
       n_mid=n+nodes_extra
       n_temp_end=n+(nodes_extra*2)
    
        !increase y nodes---------------------------- 
        do i=1,n
            do j=1,m_temp_start
                TEMP_2(1,i,j)=Nonvalue
                TEMP_2(2,i,j)=Nonvalue
            end do
        end do
     
        do i=1, n
          do j= m_temp_start+1,m_mid
              TEMP_2(1,i,j)=BOUNDING(1,i,(j-nodes_extra))
              TEMP_2(2,i,j)=BOUNDING(2,i,(j-nodes_extra))
          end do
        end do
        
        do i=1,n
             do j=m_mid+1,m_temp_end 
                TEMP_2(1,i,j)=nonvalue
                TEMP_2(2,i,j)=nonvalue
            end do
        end do
       !--------------------------------------------
       
       !increase x nodes---------------------------- 
        do i=1,n_temp_start
            do j=1, m_temp_end
                TEMP_2(3,i,j)=nonvalue
                TEMP_2(4,i,j)=nonvalue
            end do
        end do
        
        do i=n_temp_start+1, n_mid
          do j= 1,m_temp_end 
              TEMP_2(3,i,j)=TEMP_2(1,(i-nodes_extra),j)
              TEMP_2(4,i,j)=TEMP_2(2,(i-nodes_extra),j)
          end do
        end do
        
        do i=n_mid+1,n_temp_end 
            do j=1,m_temp_end 
                TEMP_2(3,i,j)=nonvalue
                TEMP_2(4,i,j)=nonvalue
            end do
        end do      
        !------------------------------------------
        n=n_temp_end
        m=m_temp_end
    
        do i=1,n
            do j=1,m
                BOUNDING(1,i,j)=TEMP_2(3,i,j)
                BOUNDING(2,i,j)=TEMP_2(4,i,j)
            end do
        end do
     
end if 

    deallocate (TEMP_2)
    
return
end subroutine increase_grid_area
!-------------------------------------------------------------------------------------------------------------------------------
!corresponding subroutine that increases the resolution of the grid on which the input surface are modelled
subroutine increase_nodes(n, m, BOUNDING, nonvalue, n_regrid, m_regrid, NODE_CHANGE, numsurf, MORE_NODES_BOUNDING, NODE_CHANGE_1)
    
    interface
     subroutine more_nodes (BIGbis, nbis, mbis, n_regridbis, m_regridbis, NODE_CHANGEbis,Nonvalue, numsurf)
     integer :: n_bis,m_bis, n_regridbis, m_regridbis, numsurf
     real :: Nonvalue
     real, allocatable :: BIGbis(:,:)
     real, allocatable :: NODE_CHANGEbis(:,:)
     end subroutine
    end interface

    integer :: n, m, i, j, numsurf, n_regrid, m_regrid, n_regrid_1, m_regrid_1
    real :: nonvalue
    real, allocatable :: BOUNDING(:,:,:)
    real, allocatable :: MORE_NODES_BOUNDING(:,:)
    real, allocatable :: NODE_CHANGE (:,:)      
    real, allocatable :: NODE_CHANGE_1(:,:)
    
     do i=1,n
        do j=1,m
           MORE_NODES_BOUNDING(i,j)=BOUNDING(numsurf,i,j)
        end do
     end do
     
     call more_nodes (MORE_NODES_BOUNDING, n, m, n_regrid_1, m_regrid_1, NODE_CHANGE_1,nonvalue,numsurf)
     call more_nodes (NODE_CHANGE_1, n_regrid_1, m_regrid_1, n_regrid, m_regrid, NODE_CHANGE,nonvalue,numsurf)
                          
return
end subroutine increase_nodes
!-------------------------------------------------------------------------------------------------------------------------------
!corresponding subroutine that rotates the user defined bounfding surfaces to match the user defined
!progradation direction of clinoform surfaces (Chapter 2, section 2.3.4) 
subroutine rotate_surfaces(n_regrid, m_regrid, NODE_CHANGE, numsurf, ROTATION, xmin, xmax, ymin, ymax, user_angle, nonvalue, resolx)

    interface
     subroutine triangulation (PointM, Nonvalue,NODE_CHANGE, n_regrid, m_regrid, resolx, resoly, xmin, xmax, ymin, ymax)
     integer :: n_regrid, m_regrid
     real :: Nonvalue, xmin, xmax, ymin, ymax, user_angle, resolx, resoly
     real :: PointM(3)
     real, allocatable :: NODE_CHANGE(:,:)
     end subroutine
    end interface
    
    character :: temp_name*1000
    integer :: numsurf, n_regrid, m_regrid
    real ::  xmin, xmax, ymin, ymax, user_angle, nreal, mreal, resolx, resoly, centeri, centerj, xcenter, ycenter, pi
    real :: rotation_angle_rads, angle_degrees, theta, coordireal, coordjreal, distance, nonvalue  
    parameter (pi=3.14159265358979323846264338327950288419716939937510)
    real,allocatable :: NODE_CHANGE(:,:)
    real,allocatable :: ROTATION(:,:,:)
    real :: PointM(3)
     
     nreal=n_regrid
     mreal=m_regrid
     resoly = (ymax-ymin)/((m_regrid-1)*1.)
     resolx = (xmax-xmin)/((n_regrid-1)*1.)
     centeri = ((nreal-1.)/2.)+1
     xcenter=xmax-(centeri-1)*resolx
     centerj = ((mreal-1.)/2.)+1
     ycenter=ymax-(centerj-1)*resoly
     
     rotation_angle_rads=(user_angle*pi)/180
     angle_degrees=(rotation_angle_rads*180)/pi
    
        
    do i = 1,n_regrid
     do j=1,m_regrid
        if (i==centeri.and.j==centerj) then; 
            ROTATION(numsurf,i,j)=NODE_CHANGE(i,j) 
        else;
            distance=SQRT(((((i*1.)-(centeri*1.))*(resolx))**2)+((((j*1.)-(centerj*1.))*(resoly))**2))
            theta = atan((((j*1.)-(centerj*1.))*resoly)/(((i*1.)-(centeri*1.))*resolx))
 
        if(((j*1.)-(centerj*1.))*((i*1.)-(centeri*1.))>=0)then;
        if (((i*1.)-(centeri*1.))>=0) then;
             theta=theta
        else;
            theta=theta+pi
        end if
        else;
        if (((i*1.)-(centeri*1.))>=0) then;
            theta=theta
        else;
            theta=theta+pi
        end if
        end if

        if (theta<0) then;
            theta=theta+2.*pi
        end if
    
            theta=theta+((2*pi)-rotation_angle_rads)
            coordireal = distance*cos(theta)/(resolx)+centeri
            coordjreal = distance*sin(theta)/(resoly)+centerj

            PointM(1) = coordireal
            PointM(2) = coordjreal
     
        if (coordireal<1.or.coordireal>n_regrid.or.coordjreal<1.or.coordjreal>m_regrid)then;
            PointM(3)=nonvalue
        else;
!-------------------------------------------------------------------------------------------------------------------------------
!subroutine that calculates the new position of the bounding surfaces after rotation 
        call triangulation(PointM, nonvalue,NODE_CHANGE,n_regrid,m_regrid, resolx, resoly, xmin, xmax, ymin, ymax)
            ROTATION(numsurf,i,j)=PointM(3)
        end if
        if (PointM(3)-1<nonvalue) then;
            ROTATION(numsurf,i,j) = PointM(3)
        end if
      end if
     end do
    end do 
    
return
end subroutine rotate_surfaces
!-------------------------------------------------------------------------------------------------------------------------------
subroutine listing(SPACE)

 integer :: SPACE, i

 open(unit=3, file="listing.txt")
  do i=1,SPACE
   write(3,*) i
  end do
 close(unit=3)

 return
end subroutine listing
!-------------------------------------------------------------------------------------------------------------------------------
!creates clinoform surfaces
!first part of subroutine initiates arrays used in subsequent subroutines
subroutine create_clinoform_surfaces(ROTATION, nonvalue, n_regrid, m_regrid, resolx, xmin, xmax, ymin, ymax, user_angle, A, filepath, proportion, Petrel)

    interface
     subroutine rotate_back(CLINOFORM, AGAIN, resolx, nonvalue, n_regrid, m_regrid, xmin, xmax, ymin, ymax, user_angle, B, Petrel, filepath)
     character :: B*5, filepath*1000
     integer :: n_regrid, m_regrid, Petrel
     real :: resolx, nonvalue, xmin, xmax, ymin, ymax, user_angle
     real, allocatable :: CLINOFORM(:,:)
     real, allocatable :: AGAIN(:,:)
     end subroutine
    end interface
    
    interface
     subroutine rotate_back_bounding(ROTATION, AGAIN, resolx, nonvalue, n_regrid, m_regrid, xmin, xmax, ymin, ymax, user_angle, B, Petrel)
     character :: B*5
     integer :: n_regrid, m_regrid, Petrel
     real :: resolx, nonvalue, xmin, xmax, ymin, ymax, user_angle
     real, allocatable :: ROTATION(:,:,:)
     real, allocatable :: AGAIN(:,:)
     end subroutine
    end interface  
      
    interface
     subroutine muddraperef(passagemud,n_regrid,m_regrid,resolx,routa,routb,rina,rinb, B,nonvalue, coordicenter, coordjcenter, dimi, dimj, proportion)
     character :: B*5
     integer :: n_regrid, m_regrid, coordicenter, coordjcenter, dimi, dimj
     real :: resolx, nonvalue,routa,routb,rina,rinb, proportion
     real, allocatable :: passagemud(:,:)
     end subroutine
    end interface
    
    interface
     subroutine mud_surface(passagemud,ROTATION,ROTATION_ADJUSTED,CLINOFORM_ADJUSTED,MUD, n_regrid,m_regrid, nonvalue, h, j0, coordicenter, coordjcenter,xmin,xmax,ymin,ymax, dimi, dimj, resolx,routa,routb,rina,rinb)
     character :: temp_name*1000
     integer :: n_regrid, m_regrid, nbcontrol, h, j0, coordicenter, coordjcenter, i, j, k, dimi, dimj 
     real :: nonvalue, coordiinf, coordjinf,xmin,xmax,ymin,ymax, resolx, routa,routb,rina,rinb
     real, allocatable :: passagemud(:,:)
     real, allocatable :: CLINOFORM_ADJUSTED(:,:)
     real, allocatable :: ROTATION(:,:,:)
     real, allocatable :: ROTATION_ADJUSTED(:,:,:)
     real, allocatable :: MUD(:,:)
     end subroutine
    end interface 
    
    interface
     subroutine rotate_back_mud(MUD,AGAIN, resolx, nonvalue, n_regrid, m_regrid, xmin, xmax, ymin, ymax, user_angle, B, Petrel, filepath)
     character :: B*5, filepath*1000
     integer :: n_regrid, m_regrid, Petrel
     real :: resolx, nonvalue, xmin, xmax, ymin, ymax, user_angle
     real, allocatable :: MUD(:,:)
     real, allocatable :: AGAIN(:,:)
     end subroutine
    end interface 
    
    interface
     subroutine alter_clinoform(ROTATION, ROTATION_ADJUSTED,n_regrid, m_regrid, nonvalue, alteration_value, CLINOFORM_ADJUSTED, xmin, xmax, ymin, ymax)
     integer :: n_regrid, m_regrid 
     real :: nonvalue, alteration_value, xmin, xmax, ymin, ymax                      
     real, allocatable :: ROTATION(:,:,:)
     real, allocatable :: CLINOFORM_ADJUSTED(:,:)
     real, allocatable :: ROTATION_ADJUSTED(:,:,:)
     end subroutine
    end interface
    
   interface 
     subroutine write_zmap_file(CLINOFORM_ADJUSTED_TEMP,temp_name,n,m,k,xmin,xmax,ymin,ymax)
     character :: temp_name*1000 
     integer :: n,m,k
     real :: xmin, xmax, ymin, ymax
     real, allocatable :: CLINOFORM_ADJUSTED_TEMP(:,:,:)
     end subroutine
   end interface 
     
    character :: parameters_input*1000, dummy*1000, B*5, temp_name*1000, filepath*1000
    integer :: n, m, numsurf, n_regrid, m_regrid, h, pl, j0, A, length, compteur, limit_sup, limit_inf, nblines, k
    integer :: Petrel, coordicenter, coordjcenter, dimi, dimj 
    real :: nonvalue, resolx, routa, routb, rina, rinb, dist1, dist2, ipole1, jpole1, jpole2, ipole2, ipole3, jpole3
    real :: radiusmax, radiusmin, xmin, xmax, ymin, ymax, jpole4, ipole4, heightbis, height, r, axis, user_angle, alteration_value
    real :: proportion
    real, allocatable :: ROTATION(:,:,:)
    real, allocatable :: CLINOFORM(:,:)
    real, allocatable :: CLINOFORM_ADJUSTED(:,:)
    real, allocatable :: CLINOFORM_ADJUSTED_TEMP(:,:,:)
    real, allocatable :: ROTATION_ADJUSTED(:,:,:)
    real, allocatable :: CLINOFORM_PAST(:,:)
    real, allocatable :: AGAIN(:,:)
    real, allocatable:: passagemud(:,:)
    real, allocatable :: MUD(:,:)
    
!-------------------------------------------------------------------------------------------------------------------------------
!read in of user defined parameters fro the clinoform-modelling algorithm    
    write(*,*) "Press 1 for Petrel or 2 for RMS"
    read(*,*) Petrel 
    write(*,*)"Enter the exponent for the clinoform shape function"
    read(*,*) pl
    write(*,*)"Enter the proximal boundary entry position"
    read(*,*) axis
    write(*,*)"Enter the file path (including file name) for the clinoform parameters"
    read(*,*) parameters_input
    write(*,*)"Enter the value for bounding surfaces alteration"
    read(*,*) alteration_value
    write(*,*)"Enter the filepath where the surfaces will be exported to"
    read(*,*) filepath
    
    allocate (CLINOFORM(n_regrid,m_regrid))
    allocate (CLINOFORM_ADJUSTED(n_regrid,m_regrid))
    allocate (CLINOFORM_PAST(n_regrid,m_regrid))
    allocate (AGAIN(n_regrid,m_regrid))
    allocate (MUD(n_regrid,m_regrid))
    
   routa_past=0
   routb_past=0
   rina_past=0
   rinb_past=0
        
    open(6, file=parameters_input)
      read(6,*)dummy
    read(6,*) routa, routb, rina, rinb
    close(unit=6)
    routa= routa/resolx
    rina= rina/resolx
    routb= routb/resolx
    rinb= rinb/resolx  
    
    
     numsurf=1
     j0 = (m_regrid*axis)
    if (routa>routb.or.routb==routa) then;
     length=routa
    else if(routb>routa)then;
     length=routb
    end if
     h=1-length

     A=0
     nblines=0 
     
     !Begin clinform calc 
do while (h<n_regrid+1)
     
     !print*, h
     open(unit=6,file=parameters_input)
     read(6,*) dummy
     do i = 1,A+1
      read(6,*) routa, routb, rina, rinb
     end do
     close(unit=6)
     routa= routa/resolx
     rina= rina/resolx
     routb= routb/resolx
     rinb= rinb/resolx 
 

    if (routa>routb.or.routb==routa) then;
    length=routa
    else if(routb>routa)then;
    length=routb
    end if
   
   if(routb>routa)then;
    ipole1 = SQRT((rinb**2)-(rina**2))
    ipole1 = h +ipole1
    jpole1 = j0
    ipole2 = SQRT((rinb**2)-(rina**2))
    ipole2 = h-ipole2
    jpole2 = j0
    
    ipole3 = SQRT((routb**2)-(routa**2))
    ipole3 = h+ipole3
    jpole3 = j0
    ipole4 = SQRT((routb**2)-(routa**2))
    ipole4 = h-ipole4
    jpole4 = j0
   end if
   
   if(routa>routb.or.routb==routa)then;
    ipole1 = h
    jpole1 = SQRT((rina)**2-(rinb)**2)
    jpole1 = j0+jpole1
    ipole2 = h
    jpole2 = SQRT((rina)**2-(rinb)**2)
    jpole2 = j0-jpole2
    
    ipole3 = h
    jpole3 = SQRT((routa)**2-(routb)**2)
    jpole3 = j0+jpole3
    ipole4 = h
    jpole4 = SQRT((routa)**2-(routb)**2)
    jpole4 = j0-jpole4
   end if
    
    compteur=0
	limit_sup=h+(length)
    limit_inf=h
	if(limit_inf<1) then;
    limit_inf=1
    end if
	if(limit_sup>n_regrid.or.limit_sup<1) then;
    limit_sup=n_regrid
    end if
    
!Calculates clinform surface from 1 to limit_inf-1        
    do i= 1, limit_inf-1
      do j=1,m_regrid    
        if (ROTATION(1,i,j)>nonvalue.and.ROTATION(2,i,j)>nonvalue) then;
            dist1 = SQRT((ipole3*1.-h*1.)**2+(jpole3*1.-j*1.)**2)+SQRT((ipole4*1.-h*1.)**2+(jpole4*1.-j*1.)**2)
            dist2 = SQRT((ipole1*1.-h*1.)**2+(jpole1*1.-j*1.)**2)+SQRT((ipole2*1.-h*1.)**2+(jpole2*1.-j*1.)**2)
            if (routa>routb.or.routb==rina) then;
                if (dist2<=2*routa.and.dist1<=2*routa)then;
                CLINOFORM(i,j)=ROTATION(2,i,j)
                end if
                if(dist2>=2*rina.and.dist1<=2*routa) then;
                radiusmin=rina
                radiusmax=routa
                height=ROTATION(2,i,j)-ROTATION(1,i,j)    
                r=SQRT((((h-h)**2)*1.)+(((j0-j)**2)*1.)) 
                !(Height function Eq. 2.1 Chapter 2)  
                heightbis= height-((r-radiusmin)*height)/(radiusmax-radiusmin)
                !(Clinoform shape function Eq. 2.8 Chapter 2)
                CLINOFORM(i,j)=ROTATION(1,i,j)+((ABS((radiusmax-r)**(pl))/ABS((radiusmax-radiusmin)**(pl)))*heightbis)
                compteur=compteur+1
                end if 
                if (dist2>=2*routa.and.dist1>=2*routa) then;
                CLINOFORM(i,j)=ROTATION(1,i,j)
                end if
            else;
                if (dist2<=2*routb.and.dist1<=2*routb)then;
                CLINOFORM(i,j)=ROTATION(2,i,j)
                end if
                if(dist2>=2*rinb.and.dist1<=2*routb) then;   
                radiusmin=rina 
                radiusmax=routa
                height = ROTATION(2,i,j)-ROTATION(1,i,j)
                r = SQRT((((h-h)**2)*1.)+(((j0-j)**2)*1.))      
                heightbis= height-((r-radiusmin)*height)/(radiusmax-radiusmin)
                CLINOFORM(i,j)=ROTATION(1,i,j) +((ABS((radiusmax-r)**(pl))/ABS((radiusmax-radiusmin)**(pl)))*heightbis)
                compteur=compteur+1                        
                end if
                if (dist2>=2*routa.and.dist1>=2*routb) then;
                CLINOFORM(i,j)=ROTATION(1,i,j)
                end if                  
            end if
         else;
         CLINOFORM(i,j)=nonvalue
         end if      
       end do
     end do
         
!calc clinoform surface
   do i=limit_inf,limit_sup
    do j=1,m_regrid
     if (ROTATION(1,i,j)>nonvalue.and.ROTATION(2,i,j)>nonvalue) then;   
        dist1 = SQRT((ipole3*1.-i*1.)**2+(jpole3*1.-j*1.)**2)+SQRT((ipole4*1.-i*1.)**2+(jpole4*1.-j*1.)**2)
        dist2 = SQRT((ipole1*1.-i*1.)**2+(jpole1*1.-j*1.)**2)+SQRT((ipole2*1.-i*1.)**2+(jpole2*1.-j*1.)**2)
        if (routa>routb.or.routb==routa) then;
           if (dist2<=2*routa.and.dist1<=2*routa)then;
           CLINOFORM(i,j)=ROTATION(2,i,j)
           end if
           if(dist2>=2*rina.and.dist1<=2*routa) then;
           !print*, dist1, dist2
           radiusmin=(rina*rinb)/(SQRT(((rina**2)*(((i-h)**2)*1.)+(rinb**2)*(((j-j0)**2)*1.))/(((i-h)**2+(j-j0)**2))*1.))
           radiusmax=(routa*routb)/(SQRT(((routa**2)*(((i-h)**2)*1.)+(routb**2)*(((j-j0)**2)*1.))/(((i-h)**2+(j-j0)**2))*1.))       
           height = ROTATION(2,i,j)-ROTATION(1,i,j)
           !(radius of clinofrom relative to clinoform origin Eq. 2.4 Chapter 2)
           r = SQRT((((h-i)**2)*1.)+(((j0-j)**2)*1.)) 
           !(Height function Eq. 2.1 Chapter 2)  
           heightbis= height-((r-radiusmin)*height)/(radiusmax-radiusmin)
           !(Clinoform shape function Eq. 2.8 Chapter 2)
           CLINOFORM(i,j)=ROTATION(1,i,j)+((ABS((radiusmax-r)**(pl))/ABS((radiusmax-radiusmin)**(pl)))*heightbis)
           compteur=compteur+1     
           end if   
           if (dist2>=2*routa.and.dist1>=2*routa) then;
           CLINOFORM(i,j)=ROTATION(1,i,j)
           end if           
        else;     
           if (dist2<=2*routb.and.dist1<=2*routb)then;
           CLINOFORM(i,j)=ROTATION(2,i,j)
           end if
           if(dist2>=2*rinb.and.dist1<=2*routb) then; 
           !(Eq. 2.5 Chapter 2, maximum radius of top ellipse)  
           radiusmin=(rina*rinb)/(SQRT(((rina**2)*(((i-h)**2)*1.)+(rinb**2)*(((j-j0)**2)*1.))/(((i-h)**2+(j-j0)**2))*1.)) 
           !(Eq. 2.6 Chapter 2, maximum radius of base ellipse)  
           radiusmax=(routa*routb)/(SQRT(((routa**2)*(((i-h)**2)*1.)+(routb**2)*(((j-j0)**2)*1.))/(((i-h)**2+(j-j0)**2))*1.))
           height = ROTATION(2,i,j)-ROTATION(1,i,j)
           !(radius of clinofrom relative to clinoform origin Eq. 2.4 Chapter 2)
           r = SQRT((((h-i)**2)*1.)+(((j0-j)**2)*1.))
           !(Height function Eq. 2.1 Chapter 2)        
           heightbis= height-((r-radiusmin)*height)/(radiusmax-radiusmin)
           !(Clinoform shape function Eq. 2.8 Chapter 2)
           CLINOFORM(i,j)=ROTATION(1,i,j) +((ABS((radiusmax-r)**(pl))/ABS((radiusmax-radiusmin)**(pl)))*heightbis)           
           compteur=compteur+1                    
           end if
           if (dist2>=2*routa.and.dist1>=2*routb) then;
           CLINOFORM(i,j)=ROTATION(1,i,j)
           end if                  
         end if
       else;
       CLINOFORM(i,j)=nonvalue
       end if
      end do    
      if (compteur>1) then;
      nblines=nblines+1
      end if
      compteur=0
     end do
     
!Ensures all clinoform values outside volume are nonvalue
    if (limit_sup<n_regrid) then;
     do i=limit_sup+1,n_regrid
       do j=1,m_regrid
        CLINOFORM(i,j)=ROTATION(1,i,j)
       end do
     end do
    end if 
 
   if (A>1) then;
    do i=1,n_regrid
       do j=1,m_regrid
        if (CLINOFORM(i,j)>nonvalue+1)then;
        if (CLINOFORM(i,j)>CLINOFORM_PAST(i,j))then;
        CLINOFORM(i,j)=CLINOFORM_PAST(i,j)
        end if
        end if
       end do
     end do
    end if

  do i=1,n_regrid
    do j=1,m_regrid
     CLINOFORM_PAST(i,j)=CLINOFORM(i,j)
    end do
  end do
           
    do i=1,n_regrid
      do j=1,m_regrid
         CLINOFORM_ADJUSTED(i,j)=CLINOFORM(i,j)
      end do
   end do
  
if (nblines>1) then;
    A=A+1
    
   open(unit=40, file="listing.txt")
    do i=1,A
	read(40,*) B
     end do
    close(unit=40)              
                             
     !if (A>0)then;
        !print*, "Generating clinoform surface",trim(B)

   ! end if
    
    ! subroutine that rotates clinoform surface to orientation of original bounding surfaces, but honour the user specified progradation direction
    call rotate_back(CLINOFORM, AGAIN, resolx, nonvalue, n_regrid, m_regrid, xmin, xmax, ymin, ymax, user_angle, B, Petrel, filepath)
    ! modifies the clinoform surfaces if user wishes to add grid layers above or below clinoforms
    call alter_clinoform(ROTATION,ROTATION_ADJUSTED,n_regrid, m_regrid, nonvalue,alteration_value, CLINOFORM_ADJUSTED, xmin, xmax, ymin, ymax)
    ! calculates the positions of the heterogeneity along clinoform surfaces
    call muddraperef(passagemud,n_regrid,m_regrid,resolx,routa,routb,rina,rinb, B,nonvalue, coordicenter, coordjcenter, dimi, dimj, proportion)
    ! creates a surface representing the heterogeneity along clinoforms
    call mud_surface(passagemud,ROTATION,ROTATION_ADJUSTED,CLINOFORM_ADJUSTED,MUD,n_regrid,m_regrid, nonvalue, h, j0, coordicenter, coordjcenter,xmin,xmax,ymin,ymax, dimi, dimj, resolx,routa,routb,rina,rinb)   
   ! subroutine that rotates surfaces representing the heterogeneity along clinoforms to orientation of original bounding surfaces, 
   ! but honours the user specified progradation direction
    call rotate_back_mud(MUD,AGAIN, resolx, nonvalue, n_regrid, m_regrid, xmin, xmax, ymin, ymax, user_angle, B, Petrel, filepath)  
end if

call xpath(h,z,A, parameters_input)
     z=z/resolx
     h=h+z
  
    if (A>1)then;

    call ypath(j0,q, A, parameters_input)
      q=q/resoly
      j0=j0+q
    end if  
     
     nblines=0
     numsurf=numsurf+1
     
end do

    ! subroutine that rotates bounding surfaces back to their original orientation
    call rotate_back_bounding(ROTATION, AGAIN, resolx, nonvalue, n_regrid, m_regrid, xmin, xmax, ymin, ymax, user_angle, B, Petrel)

    deallocate (CLINOFORM)
    deallocate (CLINOFORM_PAST)
    deallocate (CLINOFORM_ADJUSTED)
    deallocate (AGAIN)
    deallocate (MUD)
            
return
end subroutine create_clinoform_surfaces
!-------------------------------------------------------------------------------------------------------------------------------
! subroutine that rotates clinoform surfaces to orientation of original bounding surfaces, but honours the user specified progradation direction
subroutine rotate_back(CLINOFORM, AGAIN, resolx, nonvalue, n_regrid, m_regrid, xmin, xmax, ymin, ymax, user_angle, B, Petrel, filepath)

    interface
     subroutine triangulation (PointM, Nonvalue,NODE_CHANGE, n_regrid, m_regrid, resolx, resoly, xmin, xmax, ymin, ymax)
     integer :: n_regrid, m_regrid
     real :: Nonvalue, xmin, xmax, ymin, ymax, user_angle, resolx, resoly
     real :: PointM(3)
     real, allocatable :: NODE_CHANGE(:,:)
     end subroutine
    end interface
    
    interface  
     subroutine gridback (AGAIN, m_regrid, n_regrid, nonvalue, xmin, xmax, ymin, ymax, B, Petrel, filepath)
     character :: B*5, filepath*1000
     integer :: m_regrid, n_regrid, Petrel
     real :: nonvalue, xmin, xmax, ymin, ymax
     real, allocatable :: AGAIN(:,:)
     end subroutine
    end interface
    
    character :: B*5, filepath*1000 
    integer :: n_regrid, m_regrid, numsurf, Petrel
    real ::  xmin, xmax, ymin, ymax, user_angle, nreal, mreal, resolx, resoly, centeri, centerj, xcenter, ycenter, pi
    real :: rotation_angle_rads, angle_degrees, theta, coordireal, coordjreal, distance, nonvalue
    parameter (pi=3.14159265358979323846264338327950288419716939937510)
    real,allocatable :: AGAIN(:,:)
    real,allocatable :: CLINOFORM(:,:)
    real :: PointM(3)
    
     nreal=n_regrid
     mreal=m_regrid
     resoly = (ymax-ymin)/((m_regrid-1)*1.)   
     resolx = (xmax-xmin)/((n_regrid-1)*1.)
     centeri = ((nreal-1.)/2.)+1
     xcenter=xmax-(centeri-1)*resolx
     centerj = ((mreal-1.)/2.)+1
     ycenter=ymax-(centerj-1)*resoly    
                
     rotation_angle_rads = ((360-user_angle)*pi)/180   
     angle_degrees = (rotation_angle_rads*180)/pi     
     
 do i = 1,n_regrid
  do j=1,m_regrid
   if (i==centeri.and.j==centerj) then; 
       AGAIN(i,j)= CLINOFORM(i,j) ! the height for the center point is at the same node
     else;
       distance=SQRT(((((i*1.)-(centeri*1.))*(resolx))**2)+((((j*1.)-(centerj*1.))*(resoly))**2))
       theta = atan((((j*1.)-(centerj*1.))*resoly)/(((i*1.)-(centeri*1.))*resolx))

     if(((j*1.)-(centerj*1.))*((i*1.)-(centeri*1.))>=0)then;
      if (((i*1.)-(centeri*1.))>=0) then;
       theta=theta
      else;
       theta=theta+pi
      end if
     else;
      if (((i*1.)-(centeri*1.))>=0) then;
       theta=theta
      else;
       theta=theta+pi
      end if
     end if

     if (theta<0) then;
      theta=theta+2.*pi
     end if

     theta=theta+((2*pi)-rotation_angle_rads)

     coordireal = distance*cos(theta)/(resolx)+centeri
     coordjreal = distance*sin(theta)/(resoly)+centerj

     PointM(1) = coordireal
     PointM(2) = coordjreal


    if (coordireal<1.or.coordireal>n_regrid.or.coordjreal<1.or.coordjreal>m_regrid)then;
     PointM(3)=Nonvalue
    else;     
!subroutine that calculates the new position of the bounding surfaces after rotation 
     call triangulation(PointM, Nonvalue,CLINOFORM, n_regrid,m_regrid, resolx, resoly, xmin, xmax, ymin, ymax)
     AGAIN(i,j)=PointM(3)
    end if

    if (PointM(3)-1<Nonvalue) then;
     AGAIN(i,j) = PointM(3)
    end if

   end if
  end do
 end do
 
 call gridback (AGAIN, m_regrid, n_regrid, nonvalue, xmin, xmax, ymin, ymax, B, Petrel, filepath)   
 
return
end subroutine rotate_back
!-------------------------------------------------------------------------------------------------------------------------------
!Subroutine that dictates ypath
subroutine ypath(j0, q, A, parameters_input)

 implicit none
 integer :: j0, i, A
 real :: rina, routa, rinb, routb, h, q 
 character :: parameters_input*1000, ypath_dummy*25
  
  open(unit=4, file=parameters_input)
  read(4,*) ypath_dummy
 do i = 1, A+1
  read(4,*) routa, routb, rina, rinb, h, q
 end do
  close(unit=4)
 return
end subroutine ypath
!-----------------------------------------------------------------------------------------------------------------------------
!Subroutine to dictate xpath
Subroutine xpath(h,z,A,parameters_input)

 implicit none

 integer :: h, i, A
 real :: z, rina, routa, rinb, routb 
 character :: parameters_input*1000, xpath_dummy*25

  open(unit=5, file=parameters_input)
 read(5,*) xpath_dummy
 do i = 1, A+1
  read(5,*) routa, routb, rina, rinb, z
 end do
 close(unit=5)
 return
end subroutine xpath 
!---------------------------------------------------------------------------------------------------------------------
! calculates the positions of the heterogeneity along clinoform surfaces
subroutine muddraperef(passagemud,n_regrid,m_regrid,resolx,routb_P,routa_P,rinb_P,rina_P, B,Nonvalue, coordicenter, coordjcenter, dimi, dimj, proportion)                   
                        
real, allocatable :: passagemud(:,:)
integer :: P, max,n,m,dimi,dimj,coordicenter, coordjcenter, n_regrid, m_regrid               
real:: ipole1, jpole1, ipole2, jpole2, imin, imax, jmin, jmax, coordireal, coordjreal, dist, pi,Nonvalue
integer :: A_P,A_P1,h_P, j0_P, nbpointsdoneround, nbpointsdoneline, nbpointsell,k,nbpointsline,nbpointsround
real :: routa_P, routb_P, rina_P, rinb_P, resol, paramG, paramH, paramK, proportion, ell, ref, pos, mud, angle, h_Preal,j0_Preal
real:: radius, radiusmax, major, minor,xmin,xmax,ymin,ymax,radiusmin
parameter (pi=3.14159265)
character:: name*1000, B*5, temp_name*1000

    paramG=1.00354
    paramH=9.65442
    paramK=-4.3162
    
    proportion=proportion
    ell=150
    
    n=n_regrid
    m=m_regrid

    resol=resolx
    routa_P=routa_P*resol
    routb_P=routb_P*resol
    rina_P=rina_P*resol
    rinb_P=rinb_P*resol
           
    dimi=((routa_P+SQRT((n*resol)**2+(m*resol)**2))/resol)+1
    dimj=(routb_P*2)/resol+1
    coordicenter=(1+ceiling(routa_P/resol))
    coordjcenter=(1+ceiling(routb_P/resol))
    
    xmin=-(dimi-1)*resol/2
    xmax=-xmin
    ymin=-(dimj-1)*resol/2
    ymax=-ymin
       
    if (allocated(passagemud)) then;
    deallocate(passagemud)
    end if
    allocate(passagemud(dimi,dimj))     
    
    nbpointsround=0
    nbpointsline=0
    
    do i=1,dimi
        do j=1,dimj
        if(i<=coordicenter) then;
            if(((coordicenter-i)*resol/routa_P)**2+((coordjcenter-j)*resol/routb_P)**2<=1) then;
                if(((coordicenter-i)*resol/rina_P)**2+((coordjcenter-j)*resol/rinb_P)**2<=1) then;
                passagemud(i,j)=Nonvalue
                else;
                nbpointsround=nbpointsround+1
                passagemud(i,j)=0
                end if
            else;
            passagemud(i,j)=1
            end if
        else;
            if (SQRT(((coordjcenter-j)*1.)**2)<=routb_P/resol) then; 
                if (SQRT(((coordjcenter-j)*1.)**2)<=rinb_P/resol) then;
                passagemud(i,j)=Nonvalue
                else;
                nbpointsline=nbpointsline+1
                passagemud(i,j)=0
                end if
            else;
            passagemud(i,j)=1
            end if
        end if
        end do
    end do
    
    nbpointsdoneround=0
    nbpointsdoneline=0
    
    do while(nbpointsdoneround<nbpointsround*(proportion/100))
        
        pos=2
        ref=0
        mud=1
        
        do while(ref<mud)
        call random_number(angle)
        angle=angle*pi
        call random_number(pos)
        radiusmax=routb_P*routa_P/SQRT((routb_P*sin(angle))**2+(routa_P*cos(angle))**2)
        radiusmin=rinb_P*rina_P/SQRT((rinb_P*sin(angle))**2+(rina_P*cos(angle))**2)
        !distance from clinoform pinch out, Eq. 3.2, Chapter 3)
        radius=(1-pos)*(radiusmax-radiusmin)+radiusmin
        coordireal=coordicenter-(radius)*sin(angle)/resol
        coordjreal=coordjcenter+(radius)*cos(angle)/resol
        !Frequency function Eq. 3.1, Chapter 3)
        ref=paramG/(1+exp(paramH*(pos)+paramK))
        call random_number(mud)
        end do
        
        nbpointsell=0
        
        do while (nbpointsell<10)
        
                nbpointsell=0   
                call random_number(major)
                major=ell*major*(radiusmax-radiusmin)/200
                call random_number(minor)
                minor=minor*major
            
                ipole1=(SQRT(major**2-minor**2))/resol
                jpole1=0
                ipole2=-ipole1
                jpole2=-jpole1
                imin=-(major/resol)-1
                imax=(major/resol)+1
                jmin=-(minor/resol)-1
                jmax=(minor/resol)+1
               
                do i=imin,imax
                    do j=jmin,jmax
                    dist=SQRT(((ipole1*1.-i*1.)*resol)**2+((jpole1*1.-j*1.)*resol)**2)+SQRT(((ipole2*1.-i*1.)*resol)**2+((jpole2*1.-j*1.)*resol)**2)
                    if (dist<2*major) then;
                    nbpointsell=nbpointsell+1
                    end if
                    end do
                end do
                
        end do
        
        
        jpole1=coordjreal+(ipole1)*sin(angle)
        jpole2=coordjreal+(ipole2)*sin(angle)
        ipole1=coordireal+(ipole1)*cos(angle)
        ipole2=coordireal+(ipole2)*cos(angle)
 
        
        imin=floor(coordireal-(major/resol))
        if(imin<1) then;
        imin=1
        end if
        imax=ceiling(coordireal+(major/resol))
        if(imax>dimi) then;
        imax=dimi
        end if
        jmin=floor(coordjreal-(major/resol))
        if(jmin<1) then;
        jmin=1
        end if
        jmax=ceiling(coordjreal+(major/resol))            
        if(jmax>dimj) then;
        jmax=dimj
        end if

        do i=imin,imax
            do j=jmin,jmax
                if (passagemud(i,j)<0.5.and.passagemud(i,j)-1>Nonvalue) then;
                dist=SQRT(((ipole1-i)*resol)**2+((jpole1-j)*resol)**2)+SQRT(((ipole2-i)*resol)**2+((jpole2-j)*resol)**2)
                    if (dist<2*major) then;
                    passagemud(i,j)=1
                        if (i<=coordicenter) then;
                        nbpointsdoneround=nbpointsdoneround+1
                        else;
                        nbpointsdoneline=nbpointsdoneline+1
                        end if
                    end if
                end if
            end do
        end do
    
    end do
    
    do while(nbpointsdoneline<nbpointsline*(proportion/100))
    
        pos=2
        ref=0
        mud=1
        angle=pi
        
        do while(ref<mud)
        call random_number(coordireal)
        coordireal=coordireal*(dimi-coordicenter)+coordicenter
        call random_number(coordjreal)
        coordjreal=coordjreal*(dimj-1)+1
        pos=1-(SQRT(((coordjcenter-coordjreal)*resol)**2)-rinb_P)/(routb_P-rinb_P)
        radiusmax=routb_P
        radiusmin=rinb_P
        radius=(1-pos)*(radiusmax-radiusmin)
        ref=paramG/(1+exp(paramH*pos+paramK))
        call random_number(mud)
        end do

        
        nbpointsell=0
        
        do while (nbpointsell<10)
        
                nbpointsell=0   
                call random_number(major)
                major=ell*major*(radiusmax-radiusmin)/200
                call random_number(minor)
                minor=minor*major
            
                ipole1=(SQRT(major**2-minor**2))/resol
                jpole1=0
                ipole2=-ipole1
                jpole2=-jpole1
                imin=-(major/resol)-1
                imax=(major/resol)+1
                jmin=-(minor/resol)-1
                jmax=(minor/resol)+1    
               
                do i=imin,imax
                    do j=jmin,jmax
                    dist=SQRT(((ipole1*1.-i*1.)*resol)**2+((jpole1*1.-j*1.)*resol)**2)+SQRT(((ipole2*1.-i*1.)*resol)**2+((jpole2*1.-j*1.)*resol)**2)
                    if (dist<2*major) then;
                    nbpointsell=nbpointsell+1
                    end if
                    end do
                end do
                
        end do
        
        
        jpole1=coordjreal+(ipole1)*sin(angle)
        jpole2=coordjreal+(ipole2)*sin(angle)
        ipole1=coordireal+(ipole1)*cos(angle)
        ipole2=coordireal+(ipole2)*cos(angle)              
 
        
        imin=floor(coordireal-(major/resol))
        if(imin<1) then;
        imin=1
        end if
        imax=ceiling(coordireal+(major/resol))
        if(imax>dimi) then;
        imax=dimi
        end if
        jmin=floor(coordjreal-(major/resol))
        if(jmin<1) then;
        jmin=1
        end if
        jmax=ceiling(coordjreal+(major/resol))            
        if(jmax>dimj) then;
        jmax=dimj
        end if
               
        do i=imin,imax
            do j=jmin,jmax
                if (passagemud(i,j)<0.5.and.passagemud(i,j)-1>Nonvalue) then;
                dist=SQRT(((ipole1-i)*resol)**2+((jpole1-j)*resol)**2)+SQRT(((ipole2-i)*resol)**2+((jpole2-j)*resol)**2)
                    if (dist<2*major) then;
                    passagemud(i,j)=1
                        if (i<=coordicenter) then;
                        nbpointsdoneround=nbpointsdoneround+1
                        else;
                        nbpointsdoneline=nbpointsdoneline+1
                        end if
                    end if
                end if
            end do
        end do
    
    end do 
    
 return 
 end subroutine muddraperef 
 !---------------------------------------------------------------------------------------------------------------------
 ! creates a surface representing the heterogeneity along clinoforms
subroutine mud_surface(passagemud,ROTATION,ROTATION_ADJUSTED,CLINOFORM_ADJUSTED,MUD, n_regrid,m_regrid, nonvalue, h, j0, coordicenter, coordjcenter,xmin,xmax,ymin,ymax, dimi, dimj, resolx,routa,routb,rina,rinb)
        
    character :: temp_name*1000
    integer :: n_regrid, m_regrid, nbcontrol, h, j0, coordicenter, coordjcenter, i, j, k, dimi, dimj, c
    integer :: icoordinate, jcoordinate,last_value, icoordinate_2, jcoordinate_2, fci, fcj, lci, lcj
    real :: nonvalue,xmin,xmax,ymin,ymax, coordiinf, coordjinf, resolx,routa,routb,rina,rinb, radiusmax, radiusmin
    real, allocatable :: passagemud(:,:)
    real, allocatable :: CLINOFORM_ADJUSTED(:,:)
    real, allocatable :: ROTATION_ADJUSTED(:,:,:)
    real, allocatable :: MUD(:,:)
    real, allocatable :: MUD_TEMP(:,:,:)
    real, allocatable :: ROTATION(:,:,:)
     
  do i=1,n_regrid
   do j=1,m_regrid
    coordiinf=floor(coordicenter*1.-(i*1.-h*1.)) 
    coordjinf=floor(coordjcenter*1.-(j*1.-j0*1.))
    if (coordiinf<1.or.coordjinf<1.or.coordiinf+1>dimi.or.coordjinf+1>dimj) then;
    MUD(i,j)=CLINOFORM_ADJUSTED(i,j)
    else;
    nbcontrol=passagemud(coordiinf,coordjinf)*100+passagemud(coordiinf,coordjinf+1)*1000+passagemud(coordiinf+1,coordjinf)*1+passagemud(coordiinf+1,coordjinf+1)*10
    if(nbcontrol==1001.or.nbcontrol==110.or.nbcontrol==1.or.nbcontrol==10.or.nbcontrol==100.or.nbcontrol==1000.or.nbcontrol==0) then;
    MUD(i,j)=nonvalue
    else;
    MUD(i,j)=CLINOFORM_ADJUSTED(i,j)
    end if
    end if
   end do
  end do
  
do i=1,n_regrid
 do j=1,m_regrid
 if (i>1.and.i<(n_regrid-2))then;
    if (MUD(i,j)>ROTATION(2,i,j))then;
        if(MUD(i+2,j)==nonvalue) then;
            MUD(i,j)=nonvalue
        end if
    end if
 end if
 end do
end do

deallocate(ROTATION_ADJUSTED)   

return
end subroutine mud_surface
!--------------------------------------------------------------------------------------------------------------------- 
! subroutine that rotates surfaces representing the heterogeneity along clinoforms to orientation of original bounding surfaces, 
! but honours the user specified progradation direction
subroutine rotate_back_mud(MUD, AGAIN, resolx, nonvalue, n_regrid, m_regrid, xmin, xmax, ymin, ymax, user_angle, B, Petrel, filepath)

    interface
     subroutine triangulation (PointM, Nonvalue,NODE_CHANGE, n_regrid, m_regrid, resolx, resoly, xmin, xmax, ymin, ymax)
     integer :: n_regrid, m_regrid
     real :: Nonvalue, xmin, xmax, ymin, ymax, user_angle, resolx, resoly
     real :: PointM(3)
     real, allocatable :: NODE_CHANGE(:,:)
     end subroutine
    end interface
    
    interface
     subroutine gridback_mud (AGAIN, m_regrid, n_regrid, nonvalue, xmin, xmax, ymin, ymax, B, Petrel, filepath)
     character :: B*5, filepath*1000
     integer :: m_regrid, n_regrid, Petrel
     real :: nonvalue, xmin, xmax, ymin, ymax
     real, allocatable :: AGAIN(:,:)
     end subroutine
    end interface
    
    character :: B*5, temp_name*1000, filepath*1000  
    integer :: n_regrid, m_regrid, numsurf, Petrel
    real ::  xmin, xmax, ymin, ymax, user_angle, nreal, mreal, resolx, resoly, centeri, centerj, xcenter, ycenter, pi
    real :: rotation_angle_rads, angle_degrees, theta, coordireal, coordjreal, distance, nonvalue
    parameter (pi=3.14159265358979323846264338327950288419716939937510)
    real,allocatable :: AGAIN(:,:)
    real,allocatable :: MUD(:,:)
    real :: PointM(3)
    
     nreal=n_regrid
     mreal=m_regrid
     resoly = (ymax-ymin)/((m_regrid-1)*1.)   
     resolx = (xmax-xmin)/((n_regrid-1)*1.)
     centeri = ((nreal-1.)/2.)+1
     xcenter=xmax-(centeri-1)*resolx
     centerj = ((mreal-1.)/2.)+1
     ycenter=ymax-(centerj-1)*resoly    
                
     rotation_angle_rads = ((360-user_angle)*pi)/180   
     angle_degrees = (rotation_angle_rads*180)/pi     
     
 do i = 1,n_regrid
  do j=1,m_regrid
   if (i==centeri.and.j==centerj) then; 
       AGAIN(i,j)= MUD(i,j) ! the height for the center point is at the same node
     else;
       distance=SQRT(((((i*1.)-(centeri*1.))*(resolx))**2)+((((j*1.)-(centerj*1.))*(resoly))**2))
       theta = atan((((j*1.)-(centerj*1.))*resoly)/(((i*1.)-(centeri*1.))*resolx))

     if(((j*1.)-(centerj*1.))*((i*1.)-(centeri*1.))>=0)then;
      if (((i*1.)-(centeri*1.))>=0) then;
       theta=theta
      else;
       theta=theta+pi
      end if
     else;
      if (((i*1.)-(centeri*1.))>=0) then;
       theta=theta
      else;
       theta=theta+pi
      end if
     end if

     if (theta<0) then;
      theta=theta+2.*pi
     end if

     theta=theta+((2*pi)-rotation_angle_rads)

     coordireal = distance*cos(theta)/(resolx)+centeri
     coordjreal = distance*sin(theta)/(resoly)+centerj

     PointM(1) = coordireal
     PointM(2) = coordjreal


    if (coordireal<1.or.coordireal>n_regrid.or.coordjreal<1.or.coordjreal>m_regrid)then;
     PointM(3)=Nonvalue
    else;     
     call triangulation(PointM, Nonvalue,MUD, n_regrid,m_regrid, resolx, resoly, xmin, xmax, ymin, ymax)
     AGAIN(i,j)=PointM(3)
    end if

    if (PointM(3)-1<Nonvalue) then;
     AGAIN(i,j) = PointM(3)
    end if

   end if
  end do
 end do

 call gridback_mud (AGAIN, m_regrid, n_regrid, nonvalue, xmin, xmax, ymin, ymax, B, Petrel, filepath)   
 
return
end subroutine rotate_back_mud
!----------------------------------------------------------------------------------------------------------------------------------------------------
! modifies the clinoform surfaces if user wishes to add grid layers above or below clinoforms
subroutine alter_clinoform(ROTATION,ROTATION_ADJUSTED,n_regrid, m_regrid, nonvalue, alteration_value, CLINOFORM_ADJUSTED, xmin, xmax, ymin, ymax)
     
     character :: temp_name*1000  
     integer :: n_regrid, m_regrid 
     real :: nonvalue, alteration_value, xmin, xmax, ymin, ymax                      
     real, allocatable :: ROTATION(:,:,:)
     real, allocatable :: CLINOFORM_ADJUSTED(:,:)
     real, allocatable :: CLINOFORM_ADJUSTED_TEMP(:,:,:)
     real, allocatable :: ROTATION_ADJUSTED(:,:,:)

allocate(ROTATION_ADJUSTED(2,n_regrid,m_regrid)) 

do i=1, n_regrid
 do j=1,m_regrid
    if(ROTATION(1,i,j)>nonvalue)then;
        ROTATION_ADJUSTED(1,i,j)=ROTATION(1,i,j)+alteration_value
        ROTATION_ADJUSTED(2,i,j)=ROTATION(2,i,j)-alteration_value
    end if
    if(ROTATION(1,i,j)==nonvalue)then;
        ROTATION_ADJUSTED(1,i,j)=nonvalue
        ROTATION_ADJUSTED(2,i,j)=nonvalue
    end if
 end do
end do

allocate(CLINOFORM_ADJUSTED_TEMP(1,n_regrid,m_regrid))

do i=1, n_regrid
 do j=1,m_regrid
    if (CLINOFORM_ADJUSTED(i,j)==ROTATION(2,i,j)) then;
        CLINOFORM_ADJUSTED(i,j)=ROTATION_ADJUSTED(2,i,j)
    end if
    if (CLINOFORM_ADJUSTED(i,j)==ROTATION(1,i,j)) then;
        CLINOFORM_ADJUSTED(i,j)=ROTATION_ADJUSTED(1,i,j)
    end if
 end do
end do

do i=1,n_regrid
 do j=1,m_regrid
    CLINOFORM_ADJUSTED_TEMP(1,i,j)=CLINOFORM_ADJUSTED(i,j)
 end do
end do


do i=1, (n_regrid-6)
 do j=1,m_regrid
 if (i>1.and.i<n_regrid)then;
    if (CLINOFORM_ADJUSTED(i,j)==ROTATION_ADJUSTED(2,i,j)) then;
        if(CLINOFORM_ADJUSTED(i-1,j)==ROTATION_ADJUSTED(2,i-1,j).and.CLINOFORM_ADJUSTED(i+1,j)==ROTATION_ADJUSTED(2,i+1,j)) then;
            CLINOFORM_ADJUSTED_TEMP(1,i,j)=nonvalue
        end if
     end if
     if (CLINOFORM_ADJUSTED(i,j)==ROTATION_ADJUSTED(1,i,j)) then;
            CLINOFORM_ADJUSTED_TEMP(1,i+6,j)=nonvalue
    end if
 end if
 end do
end do

do i=1,n_regrid
 do j=1,m_regrid
    CLINOFORM_ADJUSTED(i,j)=CLINOFORM_ADJUSTED_TEMP(1,i,j)
 end do
end do

!temp_name="clino_adjusted.zmap+"
!k=1
!call write_zmap_file(CLINOFORM_ADJUSTED_TEMP,temp_name,n_regrid,m_regrid,k,xmin,xmax,ymin,ymax)

deallocate(CLINOFORM_ADJUSTED_TEMP)

return
end subroutine
!-------------------------------------------------------------------------------------------------------------------------------
!subroutine to rotate bounding surfaces back to their original configuration
subroutine rotate_back_bounding(ROTATION, AGAIN, resolx, nonvalue, n_regrid, m_regrid, xmin, xmax, ymin, ymax, user_angle, B, Petrel)

    interface
     subroutine triangulation (PointM, Nonvalue,NODE_CHANGE, n_regrid, m_regrid, resolx, resoly, xmin, xmax, ymin, ymax)
     integer :: n_regrid, m_regrid
     real :: Nonvalue, xmin, xmax, ymin, ymax, user_angle, resolx, resoly
     real :: PointM(3)
     real, allocatable :: NODE_CHANGE(:,:)
     end subroutine
    end interface
    
    interface
     subroutine gridback_bounding (AGAIN, m_regrid, n_regrid, nonvalue, xmin, xmax, ymin, ymax, k, Petrel)
     character :: B*5
     integer :: m_regrid, n_regrid, Petrel, k
     real :: nonvalue, xmin, xmax, ymin, ymax
     real, allocatable :: AGAIN(:,:)
     end subroutine
    end interface
    
    character :: B*5 
    integer :: n_regrid, m_regrid, numsurf, Petrel,k
    real ::  xmin, xmax, ymin, ymax, user_angle, nreal, mreal, resolx, resoly, centeri, centerj, xcenter, ycenter, pi
    real :: rotation_angle_rads, angle_degrees, theta, coordireal, coordjreal, distance, nonvalue
    parameter (pi=3.14159265358979323846264338327950288419716939937510)
    real,allocatable :: AGAIN(:,:)
    real,allocatable :: TEMP_BOUND(:,:)
    real,allocatable :: ROTATION(:,:,:)
    real :: PointM(3)
    
     nreal=n_regrid
     mreal=m_regrid
     resoly = (ymax-ymin)/((m_regrid-1)*1.)   
     resolx = (xmax-xmin)/((n_regrid-1)*1.)
     centeri = ((nreal-1.)/2.)+1
     xcenter=xmax-(centeri-1)*resolx
     centerj = ((mreal-1.)/2.)+1
     ycenter=ymax-(centerj-1)*resoly    
                
     rotation_angle_rads = ((360-user_angle)*pi)/180   
     angle_degrees = (rotation_angle_rads*180)/pi     
 
 allocate(TEMP_BOUND(n_regrid,m_regrid))
 
 do k=1,2
   do i=1,n_regrid
      do j=1,m_regrid
      TEMP_BOUND(i,j)=ROTATION(k,i,j)
      end do
   end do
     
 do i = 1,n_regrid
  do j=1,m_regrid
   if (i==centeri.and.j==centerj) then; 
       AGAIN(i,j)=TEMP_BOUND(i,j) ! the height for the center point is at the same node
     else;
       distance=SQRT(((((i*1.)-(centeri*1.))*(resolx))**2)+((((j*1.)-(centerj*1.))*(resoly))**2))
       theta = atan((((j*1.)-(centerj*1.))*resoly)/(((i*1.)-(centeri*1.))*resolx))

     if(((j*1.)-(centerj*1.))*((i*1.)-(centeri*1.))>=0)then;
      if (((i*1.)-(centeri*1.))>=0) then;
       theta=theta
      else;
       theta=theta+pi
      end if
     else;
      if (((i*1.)-(centeri*1.))>=0) then;
       theta=theta
      else;
       theta=theta+pi
      end if
     end if

     if (theta<0) then;
      theta=theta+2.*pi
     end if

     theta=theta+((2*pi)-rotation_angle_rads)

     coordireal = distance*cos(theta)/(resolx)+centeri
     coordjreal = distance*sin(theta)/(resoly)+centerj

     PointM(1) = coordireal
     PointM(2) = coordjreal


    if (coordireal<1.or.coordireal>n_regrid.or.coordjreal<1.or.coordjreal>m_regrid)then;
     PointM(3)=Nonvalue
    else;     
     call triangulation(PointM, Nonvalue,TEMP_BOUND, n_regrid,m_regrid, resolx, resoly, xmin, xmax, ymin, ymax)
     AGAIN(i,j)=PointM(3)
    end if

    if (PointM(3)-1<Nonvalue) then;
     AGAIN(i,j) = PointM(3)
    end if

   end if
  end do
 end do
 
 call gridback_bounding(AGAIN, m_regrid, n_regrid, nonvalue, xmin, xmax, ymin, ymax, k, Petrel) 
 
end do  
 
return
end subroutine rotate_back_bounding
!-----------------------------------------------------------------------------------------------------------------------------------------------------------------
!alters the clinoform surface to reflect post depositional erosion or truncation
subroutine post_processing(A, processing, base_name, top_name, model, filepath, Petrel)

interface
     subroutine new_list(A)
     integer :: A
     end subroutine
end interface

interface
     subroutine clinoform_filename(clinoform_filename_S, S, A, B, filepath)
     character :: clinoform_filename_S*1000, B*5, filepath*1000
     integer :: A, S
     end subroutine
end interface

interface
     subroutine mud_filename(mud_filename_S, S, A, B, filepath)
     character :: mud_filename_S*1000, B*5, filepath*1000
     integer :: A, S
     end subroutine
end interface

interface
     subroutine bounding_again(BOUNDING,base_name,top_name, n, m, nonvalue, xmin, xmax, ymin, ymax)
     character :: base_name*1000, top_name*1000
     integer :: n,m
     real :: nonvalue, xmin, xmax, ymin, ymax
     real,allocatable :: BOUNDING(:,:,:)
     end subroutine
end interface

interface
     subroutine increase_grid_area(n,m,BOUNDING, nonvalue, xmin, xmax, ymin, ymax)
     integer :: n,m
     real :: xmin, xmax, ymin, ymax, nonvalue
     real, allocatable :: BOUNDING(:,:,:)
     end subroutine
end interface

interface
     subroutine erosional_bounding(EROSIONAL, n, m, processing, nonvalue, xmin, xmax, ymin, ymax)
     integer :: n,m, processing
     real :: nonvalue, xmin, xmax, ymin, ymax
     real,allocatable :: EROSIONAL(:,:,:)
     end subroutine
end interface

interface
     subroutine clinoform_surf(clinoform_filename_S, S, A, B, CLINOFORM, n, m,xmin,xmax,ymin,ymax,nonvalue) 
     character :: clinoform_filename_S*1000, B*5
     integer :: n,m,A, S
     real :: xmin,xmax,ymin,ymax,nonvalue
     real,allocatable :: CLINOFORM(:,:,:)
     end subroutine
end interface

interface
     subroutine modify_surfaces(EROSIONAL,CLINOFORM_NEW, n, m, processing, clinoform_filename_S,xmin,xmax,ymin,ymax,nonvalue, Petrel) 
     character :: clinoform_filename_S*1000
     integer :: n,m, processing, Petrel
     real :: xmin,xmax,ymin,ymax, nonvalue
     real,allocatable :: EROSIONAL(:,:,:)
     real,allocatable :: CLINOFORM_NEW(:,:,:)
     end subroutine
end interface

interface
     subroutine new_bounding(EXTENT, n, m, nonvalue, xmin, xmax, ymin, ymax, model)
     integer :: n,m, model
     real :: nonvalue, xmin, xmax, ymin, ymax
     real,allocatable :: EXTENT(:,:,:)
     end subroutine
end interface

interface
     subroutine model_extent_surface(EXTENT,CLINOFORM_NEW, n, m, model, nonvalue,clinoform_filename_S,xmin,xmax,ymin,ymax)
     character :: clinoform_filename_S*1000
     integer :: n,m, model
     real :: nonvalue
     real :: xmin,xmax,ymin,ymax
     real,allocatable :: EXTENT(:,:,:)
     real,allocatable :: CLINOFORM_NEW(:,:,:)
     end subroutine
end interface

interface
     subroutine final_qc(BOUNDING,CLINOFORM,CLINOFORM_NEW,n, m,clinoform_filename_S,xmin,xmax,ymin,ymax,nonvalue, Petrel) 
     character :: clinoform_filename_S*1000
     integer :: n,m,Petrel
     real :: xmin,xmax,ymin,ymax,nonvalue
     real,allocatable :: BOUNDING(:,:,:)
     real,allocatable :: CLINOFORM(:,:,:)
     real,allocatable :: CLINOFORM_NEW(:,:,:)
     end subroutine
end interface

character :: clinoform_filename_S*1000, B*5,base_name*1000, top_name*1000,mud_filename_S*1000, filepath*1000, temp_name*1000
integer :: A, S, processing, n, m, model, Petrel
real :: xmin,xmax,ymin,ymax,nonvalue
real,allocatable :: BOUNDING(:,:,:)
real,allocatable :: EROSIONAL(:,:,:)
real,allocatable :: EXTENT(:,:,:)
real,allocatable :: CLINOFORM(:,:,:)
real,allocatable :: CLINOFORM_NEW(:,:,:)
real,allocatable :: MUD(:,:,:)
real,allocatable :: MUD_NEW(:,:,:)


call bounding_again(BOUNDING, base_name, top_name, n, m, nonvalue, xmin, xmax, ymin, ymax)
call increase_grid_area(n, m, BOUNDING, nonvalue, xmin, xmax, ymin, ymax)
if(processing==1.or.processing==2.or.processing==3)then;
        call erosional_bounding(EROSIONAL, n, m, processing, nonvalue, xmin, xmax, ymin, ymax)
        call increase_grid_area(n, m, EROSIONAL, nonvalue, xmin, xmax, ymin, ymax)
end if
if(model==1.or.model==2)then;
        call new_bounding(EXTENT, n, m, nonvalue, xmin, xmax, ymin, ymax, model)
        !call increase_grid_area(n, m, EXTENT, nonvalue, xmin, xmax, ymin, ymax)
end if
call new_list(A)

allocate(CLINOFORM(1,n,m))
allocate(CLINOFORM_NEW(1,n,m))

do S=1,A
        call clinoform_filename(clinoform_filename_S, S, A, B, filepath)
        call clinoform_surf(clinoform_filename_S, S, A, B, CLINOFORM, n, m,xmin,xmax,ymin,ymax,nonvalue)
     do i=1,n
        do j=1,m
           CLINOFORM_NEW(1,i,j)=CLINOFORM(1,i,j)
         end do
     end do  
    if(processing==1.or.processing==2.or.processing==3)then;
        call modify_surfaces(EROSIONAL,CLINOFORM_NEW, n, m, processing,clinoform_filename_S,xmin,xmax,ymin,ymax,nonvalue, Petrel)
    end if
          k=1
     temp_name=clinoform_filename_S
     !call write_zmap_file(CLINOFORM_NEW,temp_name,n,m,k,xmin,xmax,ymin,ymax)
        call final_qc(BOUNDING,CLINOFORM,CLINOFORM_NEW,n,m,clinoform_filename_S,xmin,xmax,ymin,ymax,nonvalue, Petrel)
    if(model==1.or.model==2)then;
        call model_extent_surface(EXTENT,CLINOFORM_NEW, n, m, model, nonvalue,clinoform_filename_S,xmin,xmax,ymin,ymax)
    end if
end do

deallocate(CLINOFORM)
deallocate(CLINOFORM_NEW)

allocate(MUD(1,n,m))
allocate(MUD_NEW(1,n,m))

do S=1,A
        call mud_filename(mud_filename_S,S,A,B, filepath)
        call clinoform_surf(mud_filename_S,S,A,B, MUD,n,m,xmin,xmax,ymin,ymax,nonvalue)
    do i=1,n
        do j=1,m
            MUD_NEW(1,i,j)=MUD(1,i,j)
        end do
    end do
    if(processing==1.or.processing==2.or.processing==3)then;
    clinoform_filename_S = mud_filename_S 
        call modify_surfaces(EROSIONAL,MUD_NEW, n, m, processing,clinoform_filename_S,xmin,xmax,ymin,ymax,nonvalue, Petrel)
    end if
        call final_qc(BOUNDING,MUD,MUD_NEW,n,m,mud_filename_S,xmin,xmax,ymin,ymax,nonvalue, Petrel)
end do            

if(processing==1.or.processing==2.or.processing==3)then;
        deallocate(EROSIONAL) 
 end if
if(model==1.or.model==2)then;
        deallocate(EXTENT)
end if
deallocate(BOUNDING)
deallocate(MUD)
deallocate(MUD_NEW)

close(unit=4, STATUS='DELETE')

return
end subroutine post_processing
!---------------------------------------------------------------------------------------------------------------------------------------------------------------------
subroutine new_list(A)

integer :: A, i

open(unit=4, file="Clinoform List.txt")
    do i=1, A
        write(4,*) i
    end do
close(unit=4)

return
end subroutine new_list
!----------------------------------------------------------------------------------------------------------------------------------------------------------------------
subroutine clinoform_filename(clinoform_filename_S, S, A, B, filepath)

    character :: clinoform_filename_S*1000,B*5, filepath*1000 
    integer :: i, S, A
          
    open (unit=4,file="Clinoform List.txt")
        do i=1,S
	       read(4,*) B
        end do 
    close(unit=4)
    clinoform_filename_S=trim(filepath)//"Clinoform_Surface"//trim(B)//".zmap+"
        
return
end subroutine clinoform_filename
!------------------------------------------------------------------------------------------------------------------------------------------------------------------------
subroutine mud_filename(mud_filename_S, S, A, B, filepath)

    character :: mud_filename_S*1000,B*5, filepath*1000 
    integer :: i, S, A
          
    open (unit=4,file="Clinoform List.txt")
        do i=1,S
	       read(4,*) B
        end do 
    close(unit=4)
    mud_filename_S=trim(filepath)//"Heterogeneity_Surface"//trim(B)//".zmap+"
        
return
end subroutine mud_filename
!------------------------------------------------------------------------------------------------------------------------------------------------------------------------
!re-imports bounding surfaces if post-processing is required
subroutine bounding_again(BOUNDING,base_name, top_name, n, m, nonvalue, xmin, xmax, ymin, ymax)

interface
   subroutine array_allocation (temp_name,TEMP,n,m, nonvalue,xmin, xmax, ymin, ymax, passage)
    character :: temp_name*1000 
    integer :: n,m
    real :: xmin, xmax, ymin, ymax, nonvalue
    real, allocatable :: TEMP(:,:,:)
    real, allocatable :: passage(:,:)
   end subroutine
 end interface

character :: clinoform_filename_S*1000, B*5,base_name*1000, top_name*1000, temp_name*1000
integer :: A, S, processing,n,m, SPACE
parameter (SPACE=6100)
real :: xmin, xmax, ymin, ymax, nonvalue
real, allocatable :: TEMP(:,:,:)
real, allocatable :: passage(:,:)
real,allocatable :: BOUNDING(:,:,:)

allocate (passage(SPACE,5))
allocate (TEMP(2,SPACE,SPACE))

    temp_name=base_name
    call array_allocation (temp_name,TEMP,n,m, nonvalue,xmin, xmax, ymin, ymax, passage)  
      allocate (BOUNDING(2,SPACE,SPACE))
        k=1
        do i=1,n
            do j=1,m
               BOUNDING(k,i,j)=TEMP(1,i,j)
               if (BOUNDING(k,i,j)>nonvalue.and.BOUNDING(k,i,j)<0) then;
                    BOUNDING(k,i,j)= BOUNDING(k,i,j)*-1
               end if
            end do
        end do
deallocate (TEMP)
     
allocate (TEMP(2,SPACE,SPACE))

    temp_name=top_name
    call array_allocation(temp_name,TEMP,n,m, nonvalue,xmin, xmax, ymin, ymax, passage)
    k=2
        do i=1,n
            do j=1,m
               BOUNDING(k,i,j)=TEMP(1,i,j)
               if (BOUNDING(k,i,j)>nonvalue.and.BOUNDING(k,i,j)<0) then;
                    BOUNDING(k,i,j)= BOUNDING(k,i,j)*-1
               end if
            end do
        end do
        
deallocate (TEMP)    
deallocate (passage)   

return
end subroutine bounding_again
!------------------------------------------------------------------------------------------------------------------------------------------------------------------------
!processes inputted erosinal surface into array
subroutine erosional_bounding(EROSIONAL, n, m, processing, nonvalue, xmin, xmax, ymin, ymax)

interface
   subroutine array_allocation (temp_name,TEMP,n,m, nonvalue,xmin, xmax, ymin, ymax, passage)
    character :: temp_name*1000 
    integer :: n,m
    real :: xmin, xmax, ymin, ymax, nonvalue
    real, allocatable :: TEMP(:,:,:)
    real, allocatable :: passage(:,:)
   end subroutine
 end interface

character :: erosional1_name*1000, erosional2_name*1000, temp_name*1000
integer :: processing,n,m, SPACE
parameter (SPACE=6100)
real :: xmin, xmax, ymin, ymax, nonvalue
real, allocatable :: TEMP(:,:,:)
real, allocatable :: passage(:,:)
real,allocatable :: EROSIONAL(:,:,:)

allocate (passage(SPACE,5))
allocate (TEMP(2,SPACE,SPACE))

if(processing==1.or.processing==2.or.processing==3)then;
    write(*,*)"Enter the file path (including file name) of the post-processing surface:"
    read(*,*) erosional1_name
    call array_allocation (erosional1_name,TEMP,n,m, nonvalue,xmin, xmax, ymin, ymax, passage)  
      allocate (EROSIONAL(2,SPACE,SPACE))
        k=1
        do i=1,n
            do j=1,m
                EROSIONAL(k,i,j)=TEMP(1,i,j)
                if (EROSIONAL(k,i,j)>nonvalue.and.EROSIONAL(k,i,j)<0) then;
                    EROSIONAL(k,i,j)= EROSIONAL(k,i,j)*-1
               end if
               
            end do
        end do
deallocate (TEMP)
end if

if (processing==3)then;   
allocate (TEMP(2,SPACE,SPACE))

    write(*,*)"Enter the file path (including file name) of the other post-processing surface:"
    read(*,*) erosional2_name
    call array_allocation(erosional2_name,TEMP,n,m, nonvalue,xmin, xmax, ymin, ymax, passage)
        k=2
        do i=1,n
            do j=1,m
               EROSIONAL(k,i,j)=TEMP(1,i,j)
               if (EROSIONAL(k,i,j)>nonvalue.and.EROSIONAL(k,i,j)<0) then;
                    EROSIONAL(k,i,j)= EROSIONAL(k,i,j)*-1
               end if
            end do
        end do
        
deallocate (TEMP)
endif    
deallocate (passage)   

return
end subroutine erosional_bounding
!------------------------------------------------------------------------------------------------------------------------------------------------------------------------
!modifies surface to reflect erosion
subroutine modify_surfaces(EROSIONAL,CLINOFORM_NEW, n, m, processing, clinoform_filename_S,xmin,xmax,ymin,ymax,nonvalue, Petrel) 
   
     interface
   subroutine write_zmap_file(CLINOFORM_NEW,temp_name,n,m,k,xmin,xmax,ymin,ymax)
    character :: temp_name*1000 
    integer :: n,m,k
    real :: xmin, xmax, ymin, ymax
    real,allocatable :: CLINOFORM_NEW(:,:,:)
   end subroutine
end interface
    
      character :: clinoform_filename_S*1000, temp_name*1000    
      real :: xmin, xmax, ymin, ymax, nonvalue   
     integer :: n,m, processing, k, i, j, a, Petrel
     real,allocatable :: EROSIONAL(:,:,:)
     real,allocatable :: CLINOFORM_NEW(:,:,:) 
   
     if(processing==1.or.processing==3)then;
        k=1
        do i=1,n
          do j=1,m
             if(CLINOFORM_NEW(k,i,j)<EROSIONAL(k,i,j))then;
               CLINOFORM_NEW(k,i,j)=EROSIONAL(k,i,j)
             end if
           end do
        end do         
     end if     
     
     if(processing==2) then;
        k=1
        do i=1,n
          do j=1,m
             if(CLINOFORM_NEW(k,i,j)>EROSIONAL(k,i,j))then;
               CLINOFORM_NEW(k,i,j)=EROSIONAL(k,i,j)
             end if
           end do
        end do  
     end if
     
     
     if(processing==3) then;
        k=2
        a=1
        do i=1,n
          do j=1,m
             if(CLINOFORM_NEW(a,i,j)>EROSIONAL(k,i,j))then;
               CLINOFORM_NEW(a,i,j)=EROSIONAL(k,i,j)
             end if
           end do
        end do   
     end if 
     
return
end subroutine modify_surfaces
!------------------------------------------------------------------------------------------------------------------------------------------------------------------------
!imports clinoform surfaces
subroutine clinoform_surf(clinoform_filename_S, S, A, B, CLINOFORM, n, m,xmin,xmax,ymin,ymax, nonvalue)

interface
   subroutine array_allocation (temp_name,TEMP,n,m, nonvalue,xmin, xmax, ymin, ymax, passage)
    character :: temp_name*1000 
    integer :: n,m
    real :: xmin, xmax, ymin, ymax, nonvalue
    real, allocatable :: TEMP(:,:,:)
    real, allocatable :: passage(:,:)
   end subroutine
end interface

     character :: clinoform_filename_S*1000, B*5, temp_name*1000
     integer :: n,m, A, S, SPACE
     parameter (SPACE=6100)
     real :: xmin, xmax, ymin, ymax, nonvalue
     real,allocatable :: CLINOFORM(:,:,:)
     real, allocatable :: TEMP(:,:,:)
     real, allocatable :: passage(:,:)
     
allocate (passage(SPACE,5))
allocate (TEMP(2,SPACE,SPACE))
        
    temp_name=clinoform_filename_S
    call array_allocation (temp_name,TEMP,n,m, nonvalue,xmin,xmax,ymin,ymax, passage)
     do i=1,n
        do j=1,m
           CLINOFORM(1,i,j)=TEMP(1,i,j)
           if (CLINOFORM(1,i,j)>nonvalue.and.CLINOFORM(1,i,j)<0) then;
                    CLINOFORM(1,i,j)= CLINOFORM(1,i,j)*-1
               end if
        end do
     end do 

deallocate (TEMP)    
deallocate (passage) 

return   
end subroutine clinoform_surf
!------------------------------------------------------------------------------------------------------------------------------------------------------------------------
!final check to ensure that clinoform surfaces do not exceed bounding surface volume
subroutine final_qc(BOUNDING,CLINOFORM,CLINOFORM_NEW,n,m,clinoform_filename_S,xmin,xmax,ymin,ymax,nonvalue, Petrel)

interface
   subroutine write_zmap_file(CLINOFORM_NEW,temp_name,n,m,k,xmin,xmax,ymin,ymax)
    character :: temp_name*1000 
    integer :: n,m,k
    real :: xmin, xmax, ymin, ymax
    real,allocatable :: CLINOFORM_NEW(:,:,:)
   end subroutine
end interface

interface
   subroutine array_allocation (temp_name,TEMP,n,m, nonvalue,xmin, xmax, ymin, ymax, passage)
    character :: temp_name*1000 
    integer :: n,m
    real :: xmin, xmax, ymin, ymax, nonvalue
    real, allocatable :: TEMP(:,:,:)
    real, allocatable :: passage(:,:)
   end subroutine
 end interface

     character :: clinoform_filename_S*1000, temp_name*1000
     integer :: n,m,k,i,j,SPACE, Petrel
     parameter (SPACE=6100)
     real :: xmin,xmax,ymin,ymax, nonvalue
     real,allocatable :: BOUNDING(:,:,:)
     real,allocatable :: CLINOFORM(:,:,:)
     real,allocatable :: CLINOFORM_NEW(:,:,:)
     real, allocatable :: TEMP(:,:,:)
     real, allocatable :: passage(:,:)

allocate (passage(SPACE,5))
allocate (TEMP(2,SPACE,SPACE))

temp_name="bounding1.zmap+"
call array_allocation (temp_name,TEMP,n,m, nonvalue,xmin,xmax,ymin,ymax, passage)


     k=1
     do i=1,n
        do j=1,m
           if(CLINOFORM_NEW(k,i,j)>BOUNDING(k,i,j))then;
              CLINOFORM_NEW(k,i,j)=BOUNDING(k,i,j)
           end if
           if(CLINOFORM_NEW(k,i,j)==TEMP(k,i,j))then;
              CLINOFORM_NEW(k,i,j)=BOUNDING(k,i,j)
           end if
           if(CLINOFORM(k,i,j)==nonvalue)then;
              CLINOFORM_NEW(k,i,j)=nonvalue
           end if
           if(TEMP(k,i,j)==nonvalue)then;
              CLINOFORM_NEW(k,i,j)=nonvalue
           end if
        end do
     end do

temp_name="bounding2.zmap+"
call array_allocation (temp_name,TEMP,n,m, nonvalue,xmin,xmax,ymin,ymax, passage)
     
     k=2
     do i=1,n
        do j=1,m
           if(CLINOFORM_NEW(1,i,j)<BOUNDING(k,i,j))then;
              CLINOFORM_NEW(1,i,j)=BOUNDING(k,i,j)
           end if
           if(CLINOFORM_NEW(1,i,j)==TEMP(1,i,j))then;
              CLINOFORM_NEW(1,i,j)=BOUNDING(k,i,j)
           end if
           if(CLINOFORM(1,i,j)==nonvalue)then;
              CLINOFORM_NEW(1,i,j)=nonvalue
           end if
           if(TEMP(1,i,j)==nonvalue)then;
              CLINOFORM_NEW(1,i,j)=nonvalue
           end if
        end do
     end do
        
     k=1
     temp_name=clinoform_filename_S
     if (Petrel==1) then;
       do i=1,n
        do j=1,m
            if (CLINOFORM_NEW(1,i,j)>nonvalue.and.CLINOFORM_NEW(1,i,j)>0) then;
               CLINOFORM_NEW(1,i,j)= CLINOFORM_NEW(1,i,j)*-1
            end if
        end do
       end do
     end if
     call write_zmap_file(CLINOFORM_NEW,temp_name,n,m,k,xmin,xmax,ymin,ymax)
     
     print*, "Completed QC of ",trim(clinoform_filename_S)   

deallocate (TEMP)    
deallocate (passage)
                   
return
end subroutine final_qc
!------------------------------------------------------------------------------------------------------------------------------------------------------------------------
subroutine new_bounding(EXTENT, n, m, nonvalue, xmin, xmax, ymin, ymax, model)

interface
   subroutine array_allocation (temp_name,TEMP,n,m, nonvalue,xmin, xmax, ymin, ymax, passage)
    character :: temp_name*1000 
    integer :: n,m
    real :: xmin, xmax, ymin, ymax, nonvalue
    real, allocatable :: TEMP(:,:,:)
    real, allocatable :: passage(:,:)
   end subroutine
 end interface

character :: clinoform_filename_S*1000, B*5,base_name*1000, top_name*1000, temp_name*1000
integer :: A, S, processing,n,m, SPACE, model, TEST
parameter (SPACE=6100)
parameter (TEST=3100)
real :: xmin, xmax, ymin, ymax, nonvalue
real, allocatable :: TEMP(:,:,:)
real, allocatable :: passage(:,:)
real,allocatable :: EXTENT(:,:,:)

allocate (passage(SPACE,5))
allocate (TEMP(2,SPACE,SPACE))

if (model==1.or.model==2) then;
    write(*,*)"Enter the file name of the new model extent bounding surface:"
    read(*,*) base_name
    temp_name=base_name
    call array_allocation (temp_name,TEMP,n,m, nonvalue,xmin, xmax, ymin, ymax, passage)  
      allocate (EXTENT(2,TEST,TEST))
        k=1
        do i=1,n
            do j=1,m
               EXTENT(k,i,j)=TEMP(1,i,j)
            end do
        end do
end if
     
if(model==2) then;
    write(*,*)"Enter the file name of the other new model extent bounding surface:"
    read(*,*) top_name
    temp_name=top_name
    call array_allocation(temp_name,TEMP,n,m, nonvalue,xmin, xmax, ymin, ymax, passage)
    k=2
        do i=1,n
            do j=1,m
               EXTENT(k,i,j)=TEMP(1,i,j)
            end do
        end do
end if
        
deallocate (TEMP)    
deallocate (passage)   

return
end subroutine new_bounding
!------------------------------------------------------------------------------------------------------------------------------------------------------------------------
subroutine model_extent_surface(EXTENT,CLINOFORM_NEW, n, m, model, nonvalue,clinoform_filename_S,xmin,xmax,ymin,ymax)

interface
    subroutine write_zmap_file(CLINOFORM_NEW,temp_name,n,m,k,xmin,xmax,ymin,ymax)
    character :: temp_name*1000 
    integer :: n,m,k
    real :: xmin, xmax, ymin, ymax
    real,allocatable :: CLINOFORM_NEW(:,:,:)
   end subroutine
end interface
        
     character :: clinoform_filename_S*1000 
     character :: temp_name*1000 
     integer :: n,m, model, k, i, j
     real :: nonvalue
     real :: xmin,xmax,ymin,ymax
     real,allocatable :: EXTENT(:,:,:)
     real,allocatable :: CLINOFORM_NEW(:,:,:) 
     
     if(model==1.or.model==2)then;
       k=1
        do i=1,n
          do j=1,m
             if(CLINOFORM_NEW(k,i,j)==nonvalue)then;
               CLINOFORM_NEW(k,i,j)=EXTENT(k,i,j)
             end if
           end do
        end do         
     end if
     
    if(model==2)then;
       k=2
        do i=1,n
          do j=1,m
             if(CLINOFORM_NEW(1,i,j)==nonvalue)then;
               CLINOFORM_NEW(1,i,j)=EXTENT(k,i,j)
             end if
           end do
        end do         
     end if  
           
           
     k=1
     temp_name=clinoform_filename_S
     call write_zmap_file(CLINOFORM_NEW,temp_name,n,m,k,xmin,xmax,ymin,ymax)
     
     print*, "Building Final Surface ",trim(clinoform_filename_S)
     
       
return
end subroutine model_extent_surface